fig a

library(dplyr)
library(viridis)
library(hrbrthemes)
library(ggplot2)
library(reshape2)
library(TCseq)

raw <- read.csv("clusterData_raw.csv", header = TRUE,sep=',',encoding="utf-8")
a=(raw$intron0_1+raw$intron0_2+raw$intron0_3)/3
b=(raw$intron2_1+raw$intron2_2+raw$intron2_3)/3
c=(raw$intron6_1+raw$intron6_2+raw$intron6_3)/3
d=(raw$intron12_1+raw$intron12_2+raw$intron12_3)/3
raw$a=a
raw$b=b
raw$c=c
raw$d=d
raw=as.matrix(raw)
data=apply(raw[,21:24],2,as.numeric)
rownames(data)=raw[,1]
colnames(data)=c("0hpi","2hpi","6hpi","12hpi")
wssplot <- function(data,nc=15,seed=1234)
{
	wss <- (nrow(data)-1)*sum(apply(data,2,var))
	for (i in 2:nc) {
		set.seed(seed)
		wss[i] <- sum(kmeans(data,centers = i)$withinss)
	}
	plot(1:nc,wss,type = "b",xlab = "Number of Clusters",ylab = "Within groups Sum of Squares")
}
wssplot(data)
set.seed(123)
clust_res <- timeclust(data, algo = 'cm', k = 6, standardize = TRUE)
p <- timeclustplot(clust_res, cols =2,axis.line.size = 0.6, axis.title.size = 24, axis.text.size = 26, title.size = 24, legend.title.size = 24, legend.text.size = 26)

jpeg(file="cluster1.jpg", width=8, height=4, units="in", res=500)
p[1]
dev.off()

jpeg(file="cluster2.jpg", width=8, height=4, units="in", res=500)
p[2]
dev.off()

jpeg(file="cluster3.jpg", width=8, height=4, units="in", res=500)
p[3]
dev.off()

jpeg(file="cluster4.jpg", width=8, height=4, units="in", res=500)
p[4]
dev.off()

jpeg(file="cluster5.jpg", width=8, height=4, units="in", res=500)
p[5]
dev.off()

jpeg(file="cluster6.jpg", width=8, height=4, units="in", res=500)
p[6]
dev.off()